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Direct simulation of electron transfer using ring polymer molecular dynamics: 
Comparison with semiclassical instanton theory and exact quantum methods 
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The use of ring polymer molecular dynamics (RPMD) for the direct simulation of electron transfer (ET) 
reaction dynamics is analyzed in the context of Marcus theory, semiclassical instanton theory, and exact 
quantum dynamics approaches. For both fully atomistic and system-bath representations of condensed-phase 
ET, we demonstrate that RPMD accurately predicts both ET reaction rates and mechanisms throughout the 
normal and activationless regimes of the thermodynamic driving force. Analysis of the ensemble of reactive 
RPMD trajectories reveals the solvent reorganization mechanism for ET that is anticipated in the Marcus 
rate theory, and the accuracy of the RPMD rate calculation is understood in terms of its exact description of 
statistical fluctuations and its formal connection to semiclassical instanton theory for deep-tunneling processes. 
In the inverted regime of the thermodynamic driving force, neither RPMD nor a related formulation of 
semiclassical instanton theory capture the characteristic turnover in the reaction rate; comparison with exact 
quantum dynamics simulations reveals that these methods provide inadequate quantization of the real-time 
electronic-state dynamics in the inverted regime. 



I. INTRODUCTION 

Condensed-phase electron transfer (ET) reactions are 
central to many biological and synthetic pathways for en- 
ergy conversion and catalysis .'^^1^ The development of ac- 
curate, robust, and scalable methods for the study of such 
reactions is thus a key objective in theoretical chemistry. 
Although transition state theories and rate models for 
ET have been successfully applied in complex systems,™Sl 
methods for the direct simulation and mechanistic study 
of ET dynamics in general systems remain less fully de- 
veloped. To this end, we explore the use of ring poly- 
mer molecular dynamics (RPMD) for the description of 
prototypical ET reactions between mixed-valence tran- 
sition metal ions in water, and we compare the RPMD 
approach against benchmark semiclassical and quantum 
dynamics methods. 

Fundamental theoretical challenges in the direct sim- 
ulation of ET reactions arise due to the coupling of the 
intrinsically quantum mechanical electronic transitions 
with slower, classical motions of the surrounding envi- 
ronment. Numerous semiclassical and mixed quantum- 
classical dynamics methods have been developed for 
the inves tigati on of electronically non-adiabatic reactions 
reactions ,'2"i21 i^^t existing methods do not enable mech- 
anistic studies that are independent of dividing surface 
assumptions in general systems; nor do they yield dynam- 
ical trajectories that preserve the equilibrium Boltzmann 
distribution^fi^i' and allow for the use of rare-event sam- 
pling methodologies.'22l New methods are needed to accu- 
rately describe coupled electronic and nuclear dynamics 
and to enable the efficient and robust simulation of long 
trajectories that bridge the multiple timescales of ET re- 



actions in complex systems. 

RPMD22I is an approximate quantum dynamical 
method that is based on the imaginary-t ime p ath in- 
tegral formulation of statistical mechanics .1231211 It pro- 
vides an isomorphic classical molecular dynamics model 
for the real-time evolution of a quantum mechanical sys- 
tem. Previous applications of RPMD include studies of 
molecular liquids,'2SlEl hydrogen transfer rates,^!^ and 
tunneling processes in low-dimensional systems.l^SEll A 
key feature of the RPMD method is that it yields real- 
time molecular dynamics trajectories that preserve the 
exact quantum Bo ltzma nn distribution and exhibit time- 
reversal symmetry.lSSESl These properties allow RPMD to 
be used in combination with rare-event sampling meth- 
ods for the trajectory-based analysis of quantum mechan- 
ical tu nnelin g processes in systems involving thousands of 
atoms.'3SBni'\Ye have recently extended RPMD to describe 
electronic and nuclear dynamics, including solvated elec- 
tron diffusion'*^ and non-adiabatic electron injection into 
liquid wateri— 

In the current paper, RPMD is used to directly sim- 
ulate ET dynamics in both atomistic and system-bath 
representations for mixed-valence ET in water. The cal- 
culated rates and mechanisms are analyzed in the context 
of semiclassical and exact quantum methods. A descrip- 
tion of the employed methodologies is provided in Section 
HD and Section |III| presents the details of the atomistic 
and system-bath representations. Calculation details are 
given in Section |IV[ and a discussion of the results is 
presented in Section [V] 



II. METHODS 
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Several methods are utilized to investigate ET rates 
and mechanisms, including RPMD, semiclassical instan- 
ton theory, exact quantum-mechanical dynamics, and the 



classical Marcus rate theory for ET. These methods are 
summarized below. 



A. Ring Polymer Molecular Dynamics 

The RPMD equations of moti on for a quantized elec- 
tron and N classical particles are^^ED 



v(") = u;l 



I'ql^+i) + q("-i) _ 2q("- 



-V (Q)L/cxt 



(q("),Qi,...,Q^) (1) 



and 



1 " 

v^- = - ;7xr E ^Q. ^-t (q'"' , Qi, • • ■ , Q^), (2) 



3 r. = l 



where q*^"-* and v^") are the position and velocity vec- 
tors of a^^ ring polymer bead, Qo and Vj are the posi- 
tion and velocity vectors of the j classical particle, and 
n is the number of imaginary-time ring-polymer beads. 
The intra-bead harmonic frequency is uj„/{f3h), where 
/? is the reciprocal temperature. The masses of elec- 
tron and classical particles are m^ and Mj, respectively, 
^cxt(q'"\Qi, • ■ • ,Qw) is the potential energy function 
of the system, and q'"' — q'"'. Eqs. [ij and [2| gener- 
ate a classical dynamics that we employ as a model for 
the real-time dynamics of the system.BH In the limit of 
large n, these dynamics preserve the exact Boltzmann 
distribution.!^ 

As in classical formulations of the therma l rat e 
constant ,12211111 the RPMD rate can be expressed aP^lMI 



tRPMD = lim K{t)kTST- 



(3) 



Here, fcxsT is the transition state theory (TST) approx- 
imation for the rate for a dividing surface ^(r) = ^^, 



where ^(r) is a collective variable, r = {q' 



(1) 



(n) 



Q} 



is the full position vector for the system, and 
Q = {Qi,---jQAf} denotes the set of classical parti- 
cle positions. The prefactor, K{t), is the time-dependent 
transmission coefficient that accounts for recrossing of 
trajectories through the dividing surface. An important 
feature of RPMD is that calculated rates and mechanisms 
are independent of the choice of TST dividing s urface, a s 
in exact quantum and exact classical dynamics .'22131111] 
The TST rate in Eq. [|is calculated usinJSMIlZl 



«TST 



(27r/3)-i/2(55), 



-PAF{^t 



/_lo d^e-f>^P(0 



(4) 



Here, F{^) is the free energy (FE) along ^ 

we(r)-c.))' 



(5) 



where fj. is a reference point in the reactant region, 
ancPHlol 



Sc(r) 



■ d 

E 



1 (dgr) 



ni/2 



(6) 



The scalar r, e {r} in this equation indicates either a 
ring-polymer or classical particle degree of freedom, rrii 
is the corresponding mass, and d is the total number 
of degrees of freedom in the system. In Eqs. |4] and [5] 
(...) denotes the equilibrium ensemble average 



/dr/dve-^^" (■■'-)(...) 



(7) 



and 
ble 



Here, 



denotes the average in the constrained ensem- 



/ rfr / dv e-'5^"("-'^) (. . . )(5(e(r) - C*) 



(8) 



where rrib is the fictitious Parrinello-Rahman mass,'221 
v={v(i),...,v("),Vi,...,V^},and 

a=l 

+ -E^-t(q'"^Q) (10) 

is the full potential energy function for the ring polymer. 
The transmission coefficient in Eq. [3] is obtained from 
the flux-side correlation function using 



<t) 



ioh{art)-e) 



Co/^Uo 



(11) 



where h{^) is the Heaviside function, ^o is the initial ve- 
locity of the collective variable in an RPMD trajectory 
released from the dividing surface, and r^ is the time- 
evolved position of the system along that trajectory.l^Sl 



B. Semiclassical Instanton Theory 

The "Im F" premise in semiclassical rate theory relates 
the thermal rate constant in the deep-tunneling regime 
to the analytical continuation of the partition function 
into the complex planCjI^il'SSl 



PhQ, 



Imi?, 



(12) 



where Q,. is the reactant partition function and Im R is 
the imaginary part of the analytical continuation of the 
partition function for the full system. In the steepest- 
descent limit, the "Im F" description is equivalent to the 
flux-side time correlation formulation- of semiclassical 
instanton (SCI) theory.lSS] We adapt this approach to de- 
scribe the transfer of a single quantized electron in a clas- 
sical solvent. 

The partition function for the full system in the ring- 
polymer representation can be expressed 



Qn = cfdQi„{q) 



(13) 



where c = Jlf^i 5^^ 



I3h^ 



««'-(^)"7*'°'" 



-A({q(°)};Q)/;^ 



and 



A{{q^'^^};Q) = {mUn{r) 



(15) 



is the classical action for a periodic trajectory in imag- 
inary time. The notation presented here assumes that 
the quantized electron moves in a single dimension. At 
each solvent configuration, the steepest-descent approxi- 
mation to /„(Q) is obtained by expanding j4({g("^}; Q) 
to second order about its global minimum {9*-"^}, for 
which the electron ring-polymer coordinates obey the sta- 
tionary condition 



1 " 

-E 

n ^-^ 

a=l 



d 



ag(") 



U,: 



('z<"',Q 



,(a)=0(°) 



ojUq^"+'^ 



~ia-l) _ 2q{c 



The steepest-descent approximation yields 



/n(Q) 



1 



Vdct K 



-A({g<°)};Q)/n 



-0. (16) 



(17) 



where K is the Hessian matrix given by 



K„ 



U}„ 



92 



h dq(t'^dq('''> 



A({r'};Q) 



{g(.^)} = {q(o,)} 



(18) 



and where (det K) = Y[i=i Vi is obtained from the nor- 
mal mode frequencies, {rji}. 

For a reaction with a barrier, a saddle point satisfies 
the stationary condition in Eq. |16[ and the Hessian ma- 
trix exhibits an imaginary normal-mode frequency, 771. 
By analytically continuing 771 onto the real axis, and 
by integrating out the zero-frequency normal mode that 
is associated with the cyclic permutation of the ring- 
polymer beads, we obtain the steepest-descent SCI rate— 



«sci 



dQ2:„(Q) 



(19) 



where 

2:„(Q) 



n 1\ 1/2 

2Trh 



1 
Vdet' K 



-A({?(°)};Q)/ft 



(20) 



(det' K) = Y[ |?7iP is obtained from a product that ex- 
cludes the zero-frequency mode, and 



Sn = E('?~^"^'^-9^"^)' 



(21) 



a=l 



Formal connections between path-integral s tatist ics 
and reactive tunneling have long been recognizedJ ^^ I ^^" ^ 
In particular, Althorpe and coworkers^^ have recently 
emphasized the connection between the TST limit of 
RPMD and the re versib le action work (RAW) formula- 
tion of SCI theory.'51'^3' To the extent that Eq. Il9 



IS an 



(14) harmonic approximation to the RAW SCI formulation,'^ 



tRPMD — 



h 



SCI, 



(22) 



where a = 27r(/3/i|77i|)^^, and Kq is the transmission co- 
efficient through a dividing surface that minimizes the 
recrossing of RPMD trajectories. 



C. Exact Quantum Dynamics 

We obtain numerically exact quantum dynamics for 
ET using the Quasi-Adiabatic Path Integral method 
(QUAPI).[£Sl2H! The method is applied to a redox system 
composed of two diabatic electronic states and a coordi- 
nate representing polarization of the solvent dipole field; 
the solvent coordinate is in turn linearly coupled to a 
harmonic oscillator bath. 

The Hamiltonian for the redox systenpS 



Hs 



Ps 

2m, 



Vl2{s) V22{s) 



(23) 



where s is the solvent coordinate, Ps is the conjugate 
momentum, and rrig is the effective solvent mass. Here, 
Vii(s) and V22(s) are diabatic states corresponding to 
reactant and product states for ET, and Vi2{s) is the 
electronic coupling. The Hamiltonian describing the bath 
modes and their coupling to the solvent coordinate is 



Hf 



/ 

E 

i=i 



2Mi 



^ 1 

E 



M.uj', 



Q, 



j=i 



M,c.2 



(24) 



where M,-, lo 



J , and Qj are the mass, frequency and the 
position of the j*^ bath mode, respectively, and Cj is the 
strength of the coupling between the j^^ bath mode and 
the solvent coordinate. 

The exact quantum mechanical rate constant can be 
expressed in terms of the symmetrized real-time flux-flux 
correlation function,"^ 



1 



fco = lim -^ / C'PF(t)dt, 



t' 



(25) 



where 



Cff(0 = Ti[jFe'"<'''Te-'"'-'l% 



(26) 



and Qr is the reactant partition function. Here, 
H = Hs+Hbis the fuU ET Hamihonian, J' = jr[H, 1^2] 
is the operator for the flux between the reactant and 
product electronic states, V2 — |2)(2| is the projec- 
tion operator for the product electronic state, and 
tc = t — i(3h/2 is the complex time. The propagators are 
discretized into Af time slices of length Ate, and the trace 
in Eq. [26] is expanded to yield 



CffW= / rfQo( Qo 

X 

X {s^f+2,cr^f+2\J^\s^f+l,<7^f+l) 



Sl,(Jl\J^\s2J^+2,Cr2Af^ 



n K,sfc|e'^^**/''|afe.i,sfc_i) 



X 



Yl (o-fcjSfcl 



-iHAt^/h 



<^k-l,Sk~l) 



k=2 



Qo , (27) 



where Qo represents the bath degrees of freedom, Sk is 
the solvent coordinate, and ak is the electronic state at 
complex time slice k. 

The propagators in Eq. [27] are factorized using the 
quasi-adiabatic short-time approximation— 



iHAt^/h 



-iHBAt^/2h -iHsAt^/h -iHBAt^/2h 



(28) 



Analytical integration over the bath modes then yields 



CFF(i) = ^Re [C7i(2, 2, 1, 1; Q - ^2(2, 1, 2, 1; t. 



where 



+ C3(l,l,2,2;te)-C4(l,2,l,2;te)], 

Ci{(Ti,Cr_\f+i,<7_\f+2,<72Af+2;tc) 

= / dsi • • • / dsf^ / dsf^+2 ■ ■■ ds2Af+i 
222 2 

^E--- E E ■•• E ^'(«''^;^c 

0-2=1 (TjV = 1 ''■Ar+3 = l a-2Ar+i=l 



(29) 



(30) 



and 



/,(s,^;te) - Vi2{si) Vi2{s^+2)K{s,(T;t,)I{s). (31) 

Here, S2^f+2 = si, sjs/+2 = sj^+i, and we 
have introduced the notation s = {si, . . . , S2AA+2} and 

(T = {(Tl,. . . ,0'2A/'+2}- 

In Eq. ]31[ the path- integral expression for the complex- 
time propagators of the system Hamiltonian is given by 



2AA+2 



CTfe-ljSfc-l/ 



k=N+Z 



M+1 



(ak,Sk\e '' |crfc_i,Sfe_i) 



(32) 



TABLE I. Complex times tk used to calculate the {B^k'}- 



k 



tk 



2,...,AA+1 

AA + 2 

A/" + 3, . . . , 2A/' -H 2 

2Af + 3 





(A: - l/2)Afc 

t ~ il3h/2 

{2Af + 3/2-k)At* -iph 



The matrix elements in Eq. 32 are obtained using the 
numerically exact expression 

/ I -iHsAt^/h\ \ 

{sk,c7k\e "' |sfe_i,crfe_i) = 



Mo 

E 

771 — 1 



0m(sfc, akWmisk-i,'yk-i)e-'''"^'^'^/\ (33) 



where (j)mis, a) and Em are the eigenstates and eigenen- 
ergies of Hg, respectively, and Mq is the number of eigen- 
states included in the expansion. 

The discretized form of the non-local influence func- 
tional in Eq. [3TJ which accounts for bath-induced elec- 
tronic transitions in the system, is 



X(s) = Zo exp - 



2Af+2 k 

EE 

k=l k' = l 



Bkk'Sk Sk' 



(34) 



where Iq i s the pa rtition function of the uncoupled bath 
oscillators .iSSESEl] xhe diagonal elements of {Bkk'} de- 
scribe local contributions to the bath response function 
from a particular complex time slice k along the adia- 
batic path, and the off-diagonal elements describe non- 
local contributions. For the case of linear system-bath 
coupling, the diagonal matrix elements are given by 



Bkk- ^ 



^j{tk+l — tk) 



^ MjUj^ sinh(/3wj/2) """ V 2 

^jitk+i -tk +i/3)' 



(35) 



and the off-diagonal matrix elements are given by 



/ 



Bkk'-Yl 



sm 



^j{tk+i — tk) 



^J^MjUj^ 8mhi/3uj,/2)""'\ 2 

^jitk+i - tk'+i + tk - tk' + i(i) 



X cos 



X sm 



Wj(ifc'+i — tk') 



(36) 



The complex times tk in Eqs. 35 and 36 are provided in 
Table [I] 



D. Marcus Theory for ET in a Classical Solvent 



k=2 



In the Marcus theory for ETjEESHZll electronic transi- 
tions occur at solvent geometries for which the donor and 



acceptor electronic states are isoenergetic. In the limit 
of weak electronic coupling and classical solvent motions, 
the ET rate is thus 



27r 



'MT 



l-Ki^)""=- 



,9AG* 



(37) 



where V12 is the electronic coupling matrix element, 



AG* = 



(AG° + A)^ 
4A 



(38) 



A is the solvent reorganization energy, and —AG" is the 
thermodynamic driving force for the ET reaction. The 
rate expression in Eq. |37| exhibits three distinct regimes of 
behavior as the driving force is varied relative to A. In the 
normal regime, where — AG" < A, the rate increases with 
increasing driving force. A turnover in this trend occurs 
in the activationless regime, for which —AG" « A. In the 
inverted regime, for which -AG" > A, the rate decreases 
with increasing driving force. 

In the current study, we use implementations for Mar- 
cus theory, SCI theory, and RPMD in which the sol- 
vent degrees of freedom are treated classically; the role 
of nuclear quantum effects in diminishing the deg ree of 
turnover for the ET rate in the inverted regime^ESl jg not 
considered here. 



III. SYSTEMS 

ET dynamics is studied using both all-atom and 
system-bath representations for mixed-valence transition 
metal ions in water. These representations are described 
in this section. 



A. Atomistic Representation for ET 

The atomistic representation for the ET reaction 
(Fig. [1]) is described using the potential energy functiorpS 



c/cxt(q, Q) - c/soi(Q) + ^c-soi(q, Q) 

+ f/c-M(q, Q) + c/m-soi(Q), 



(39) 



where q is the electron position and Q is the set of N 
classical solvent atom positions. Solvent-solvent inter- 
actions, C/soi(Q), are described using the simple point 
charge (SPC) modeP^for explicit, rigid water molecules. 
The remaining interactions are described below, with the 
values of the parameters provided in Table lllj 

The electron-water interactions are described using the 
pairwise pseudopotentiaP^ 



N 



Ue-soi (r) = Y. ^-soi (rk) , 












^Jt^i 
















v>S 



i j^^^ 



(c) 
IV 



P-y-FJP^-^?^^i 







FIG. 1. Snapshots of the atomistic representation for the 
ET reaction, with the donor and acceptor metal ions shown 
in yellow, the electron ring polymer in black, and the water 
molecules in red and white. Typical configurations of the 
symmetric ET system are presented with the electron ring 
polymer (a) in transition between the redox sites, (b) in the 
reactant basin, and (c) in the product basin. 



where rj, = |q — Qfe|. For cases in which the atom index 
k corresponds to a hydrogen atom. 



qne 



ru < r" t 



Tjk 



{rk) 



qHe H 



(40) 



and when k corresponds to an oxygen atom. 



UtUrk) 



qoe 



(41) 



fc=i 



47reorfe 
Electron-ion interactions are described using 

(7e.M(q) = C/e.D(|q - QdD + C/c-A(|q - QaI), (42) 

where Qd and Qa denote the respective positions of the 
donor and acceptor metal ions, which are held fixed at 



TABLE II. Parameters for the atomistic representation of 
ET. 



Parameter 



Value 





Qa 


/A 






Qd 


/A 






' cut 


/A 






' cut 


/A 






qo 


/e 






qn 


/e 






qM 


/e 




70 


/ (kca 


1/mol 


A^) 



(0,0,-3.25) 

(0,0,3.25) 

1.0 

1.1 

-0.84 

0.42 

3.0 

6392.7 



TABLE III. Values of the asymmetry parameter e considered 
in the atomistic representation and the corresponding ther- 
modynamic driving force regimes. 



Case 

II 
III 
IV 

V 
VI 
VII 



0.0 
0.1 
0.2 
0.3 
0.4 
0.6 
0.7 



ET Regime 
Symmetric 
Normal 
Normal 
Activationless 
Inverted 
Inverted 
Inverted 



a separation of 6.5 A. These interactions are described 
using Shaw- type pairwise pseudopotentials.l^SlFor the ac- 
ceptor metal ion, 



[/e-A(r) 



(qM + e) e 
(qM + e) e 



r < r. 



M 



(43) 



r > r, 



M 



where r 



|q ~ Qa|, and for the donor metal ion, 



UMr) = 



qMe 

qMC 



r < r, 



M 
cut 



47reor 



r > r. 



M 



(44) 



with r = |q — QdI- The asymmetry parameter, e, ad- 
justs the thermodynamic driving force for the ET reac- 
tion while leaving the solvent reorganization energy un- 
changed. The values of e considered in this study and 
the corresponding ET regimes are presented in Table [iTT] 
The ion-water interactions are given by 



N 



Uu-soliQ) = Y. «sol(Qfc) ■ 



uLoiiQk)) 



(45) 



fc=l 



For cases in which atom index k corresponds to a hydro- 
gen atom, 

^jk /p, ^ _ qnqM 

and when k corresponds to an oxygen atom. 



(46) 



ULoiiQk 



70 



qoqM 



IQd — Qk 



47reo|QD -Qa 



(47) 



TABLE IV. Values of the asymmetry parameter e considered 
in the system-bath representation, the corresponding ther- 
modynamic driving force regimes, and the electronic coupling 
matrix element, V12P] 



Case 


Model SBl 

e IV12I 


Model SB2 

\V^2\ 


ET Regime 


I 


0.0 


6.6860 


0.0 


2.0662 


Symmetric 


II 


0.05 


6.4837 


-0.015 


2.0916 


Normal 


111 


0.10 


6.1300 


-0.025 


2.1088 


Normal 


IV 


0.20 


5.4840 


-0.050 


2.1524 


Activationless 


V 


0.30 


4.9120 


-0.075 


2.1971 


Inverted 


VI 


0.40 


4.4040 


-0.100 


2.2427 


Inverted 



The coupling IV12I is given in units of a.u./lO^ for Model SBl 
and a.u./lO^ for Model SB2; e is in atomic units. 



The potential energy functions associated with the accep- 
[(Q), are obtained by replacing Qq with 



tor ion, U^ 



Qa in Eqs. |46] and |47] These ion-water potential en- 
ergy functions include electrostatic interactions combined 
with short-range repulsive terms that reproduce the oc- 
tahedral coordination structure of the solvated ions.'^S 



B. System- Bath Representations for ET 

The system-bath representation for the ET reaction is 
described in the position basis using the potential energy 
function 

C/cxt(g, s, Q) = C/o-M (g) + C/o-soi (g, s) + Ub (s, Q) , (48) 

where the first two terms comprise the system potential, 
and Ub is the potential energy contribution due to the 
bath. The scalar coordinates q and s are the positions of 
the electron and the solvent mode, respectively. 

The first term in the system potential energy function 
models the ion-electron interaction, 



Uc-u{q) 







otherwise. 



(49) 

This one-dimensional (ID) potential energy function con- 
sists of two Coulombic wells capped by parabolic func- 
tions to remove the singularity; it is continuous, and its 
derivative is piecewise continuous over the full range of 
q. The coefficients in Eq. |49] are provided in Appendix 
[A] and the values of e considered for the system-bath 
representation are presented in Table |IV| 

The second term in the system potential energy func- 
tion models the solvent and its interactions with the 
transferring electron. 



Uc-soiiq, s) = /xstanh {(f>q) + 



2 2 

m^uj^ s . 



(50) 



TABLE V. Parameters for the system-bath representation of 

et13 



and 



Parameter 


Model SBl 




Model SB2 


TA/k 


3.25 




2.72435 


tb/A 


-3.25 




-2.72435 


M 


0.0230725 




0.0114265 


/ 




12 




^s 




0.00228 




LOc 




0.00228 




M 




1836.0 




nis 




1836.0 




V/Muj^ 




1.0 




(f>/rA 




3.0 





^ Parameters given in atomic units, unless otherwise specified. 



The first term on the RHS of this equation describes the 
couphng of the electronic dipole of the redox system to 
the solvent dipole, and uj^ is the effective frequency of the 
solvent coordinate. 

The harmonic oscillator bath potential in Eq. |49] has 
the same form as in Eq. [24) 



/ 



C/b(.s,Q) = ^ 



j=i 




Af.? 



(51) 



The bath exhibits Ohmic spectral density with cutoff fre- 
quency Wc, 



J{uj) = rjLoe 



(52) 



where the dimensionless parameter rj determines the 
strength of coupling between the system and the bath 
modes. ^^ The continuous spectral density is discretized 
into / oscillators with frequencies^^ 



-Wc log 



J - 0.5 
/ 



and coupling constants 



/ 2riMuj^ 



1/2 



(53) 



(54) 



where j = 1 . . . f. 

In the current paper, we use two sets of parameters 
for the system-bath representation. Model SBl is con- 
structed to reproduce the energy-scales of the atomistic 
representation, and Model SB2 uses parameters that are 
numerically less demanding for the QUAPI calculations. 
The parameters for the models are given in Table |V] 

As indicated previously, the QUAPI method is imple- 
mented using a discrete representation for the diabatic 
states of the redox system (Eq. 23 1. The system rep- 



resentation in the position basis described in Eq. [48] is 
therefore transformed to the electronic diabatic basis for 
the QUAPI calculations. The resulting diagonal matrix 
elements for the system potential energy are 



1/22(5) = 0-28^ + hs + C2, 



(56) 



and the constant off-diagonal elements V12 are reported 
in Table IIVI The details of this transformation and the 
values of the coefficients in Eqs. [55] and |56] are given in 
Appendix [B] 



IV. CALCULATION DETAILS 

A. Atomistic Representation 

The atomistic system includes 430 SPC water 
molecules in a cubic simulation cell with periodic bound- 
ary conditions. The side-length of the cell is L = 23.46 A. 
All calculations are performed at a temperature of 
T ~ 300 K, and all pairwise interactions are truncated 
at a distance of rcut = L/2. Long-range electrostatics 
are treated by the force-shifting algorithm,-^ where the 
Coulombic portion of each potential is multiplied by a 
damping function 5(r), such that both the potential and 
its derivative smoothly vanish at r — rcut- Specifically, 



5(r) = 



1 



.0, 



2r 





r < rcut 


2 ' 
'"cut Tc-at 




r > rcut 



(57) 



Vii{s) — ais^ + his + ci 



(55) 



Force-shifting reduces the unphysical structuring of water 
near the the cutoff radius ,1^ and it is found to have little 
effect on the solvent environment of the redox system. 



1. RPMD 

The atomistic RPMD simulations are implemented in 
the DL_POLY molecular dynamics package.l^In all sim- 
ulations, the RPMD equations of motion are evolved us- 
ing the velocity Verlet algorithm, ^^ and the constraints 
in the rigid-body water model are implemented using the 
RATTLE algorithm.'^ The electron is quantized with 
n = 1024 ring-polymer beads. As in previous RPMD 
simulations, each timestep for the electron ring polymer 
involves separate coordinate updates due to forces aris- 
ing from the physical potential and due to exact evolu- 
tion of the purely harmonic portion of the ring-polymer 
potential. The resulting integration algorithm is time- 
reversible and symplectic. 

Several collective variables are used to characterize the 
ET reaction in the atomistic representation. The position 
of the electron is described by a ring-polymer progress 
variable, or "bead-count" coordinate, defined as 

1 " 1 

Q = l 

where h = 1.25 A~^, and the metal ions are symmet- 
rically positioned on the z-axis. We also consider the 



solvent collective variable 



N 



At/(Q) 



4^^o .=1 



qk 



qk 



Qd ~ Qk 



IQA-Qfcl 



(59) 

where qk G {qH,qo} is the charge on solvent atom 
k. This solvent collective variable, which is familiar 
from earlier simulation studies of Marcus theory?^ ^'^ de- 
scribes the energy difference between the electonic dia- 
batic states in the tight-binding approximation. 

The RPMD rate in Eq. [3] is calculated from the prod- 
uct of the TST rate and the transmission coefficient. The 
TST rate described in Eq. H is obtained from F(/b), 
the FE profile in the bead-count coordinate. This FE 
profile is calculated using umbrella sampling and the 
weighted histogram analysis method (WHAM), as de- 
scribed below. '^^ ^^ 

For each value of the asymmetry parameter e, the fol- 
lowing umbrella sampling protocol is used. The region 
/b — 0.06 — 0.94 is sampled with 22 trajectories that are 
harmonically restrained to uniformly spaced values of /b 
using a restraint force constant of 1.195 x 10^ kcal/mol. 
Likewise, the regions /b = 0.945-1.0 and /b = 0.0-0.055 
are each sampled with 11 uniformly spaced windows us- 
ing a higher force constant of 1.195 x 10^ kcal/mol. The 
regions of /b = 0.986 - 0.991 and /b = 0.009 - 0.015 are 
each sampled with 5 uniformly spaced windows using a 
force constant of 1.195 x 10^ kcal/mol. The equilibrium 
sampling trajectories are performed using path-integral 
molecular dynamics (PIMD) with a Parrinello-Rahman 
mass of 364.6 a.u., which allows for a timestep of 0.025 fs; 
this choice of mass does not affect the calculate d FE 
profile or any other equilibrium ensemble aver age .l^SEZI 
Each sampling trajectory is run for at least 50 ps, and 
thermostatting is performed during the trajectory cal- 
culations by resampling the particle velocities from the 
Maxwell-Boltzmann (MB) distribution every 1.25 ps. 



The transmission coefficient in Eq. 11 is calculated us- 
ing RPMD trajectories that are released from the divid- 
ing surface at /^. For each value of e, the dividing sur- 
face is chosen to coincide with the maximum along the 
FE profile, i^(/b)- The positions of the dividing surfaces 



ft 



for 



are set to f^ = (0.5, 0.7, 0.8, 0.96, 0.98, 0.98, 0. 

the different e-cases (I, II, . . . , VI) in Table III Between 



400 and 1200 trajectories are released for each value of 
e. Each RPMD trajectory is evolved for 40 fs with a 
timestep of 5 x 10~^ fs and with the initial velocities 
sampled from the MB distribution. Initial configurations 
for the released RPMD trajectories are selected every 
100 fs from eight long, independent PIMD sampling tra- 
jectories that are constrained to the dividing surface £^^. 
These sampling trajectories are thermostatted by resam- 
pling the velocities every 200 fs, and the constraint to 
the dividing surface is enforced using the RATTLE algo- 
rithm. 

Two-dimensional (2D) FE surfaces in the ring- 
polymer centroid coordinate and the solvent coordinate, 
F{z,AU), are used for the analysis of the ET reaction 



mechanism. For a given value of e, the 2D FE surface is 
constructed using PIMD sampling trajectories that are 
harmonically restrained in both z and AU coordinates. 
The z coordinate is sampled using 43 uniformly spaced 
windows in the region of —3.575 A to -1-3.575 A with a 
harmonic restraint force constant of 169.7 kcal/mol A^^. 
To ensure adequate sampling of ring-polymer configura- 
tions spanning both metal ions, we use four additional 
sampling trajectories that are harmonically restrained to 
z = ±2.7625 A and z = ±2.925 A with a force con- 
stant of 452.5 kcal/mol A^^. The solvent coordinate is 
sampled with 15 uniformly spaced windows in the range 
— 130 to + 150 kcal/mol using a harmonic restraint force 
constant of 0.023 (kcal/mol)~^. Each sampling trajec- 
tory is run for at least 50 ps, with velocities resampled 
from the MB distribution every 500 fs. We note that 
/b is a good progress variable for ET throughout the en- 
tire regime of the thermodynamic driving forces, whereas 
the ring-polymer centroid is not. In the ET inverted 
regime, the centroid does not fully distinguish between 
ring-polymer configurations in the reactant and product 
basins; no such difficulty is experienced in the calcula- 
tions reported here. 



2. Marcus Theory 



Marcus 
Eqs 



37 and 38 



theory rates are calculated 
The driving force, —AG", 



usmg 
is ob- 
tained from F{AU) as the difference between the free 
energies of the reactant and product minima; these 
values are reported in Table |VI[ To the extent that the 
tight-binding approximation holds, the reorganization 
energy, A, is identical for all e, and we confirm that this 
is very nearly the case in our calculations. For the case 
of symmetric ET (e = 0), the reorganization energy is 
calculated using A = 4:F{AU)\au=o and is found to be 
69.7 ±0.7 kcal/mol. 

The coupling matrix element in Eq. 37 IV12I, is calcu- 
lated as 21^12! = £'1 — E'oj where Eq and Ei are the two 
lowest eigenenergies of the electron in the potential of 
the isolated metal ions with e = 0. These eigenenergies 
are obtained with an iterative, block Lanczos scheme,!^ 
performed on a uniform grid of 64 x 64 x 64 points span- 
ning the cubic simulation cell. The iterative Lanczos 
calculation employs 200 Krylov vectors and an expo- 
nential transform parameter of /3l = 0.1. The block 
Lanczos refinement uses ten blocks of five Krylov vec- 
tors. This yields a value for the tunnel splitting of 
IT/12I = 0.0403 kcal/mol (6.43 x 10"^ a.u.), which is con- 
sistent with previous calculations. 7- This value for the 
tunnel splitting was assumed to be insensitive to pres- 
ence of solvent, as has been previously demonstrated,'^, 
and independent of the value of the asymmetry parame- 
ter e. The validity of this latter assumption is confirmed 



for the system-bath models (see Table IV) 



B. System-Bath Representation 



2. Marcus Theory 



As in the atomistic representation, the calculations 
in the system-bath representation are performed at 
T = 300 K. The harmonic bath is discretized using 
/ = 12 modes. 



1. RPMD 

RPMD rates for the system-bath models are also cal- 
culated with the electron quantized using n — 1024 ring- 
polymer beads. For each value of e, the FE profile, F{ft,), 
is obtained from umbrella sampling along the /b coor- 
dinate. For both system-bath models, SBl and SB2, 
-F(/b) is sampled with two sets of harmonically restrained 
PIMD trajectories. The region of /b = 0.06 - 0.94 is 
sampled with 45 trajectories that are harmonically re- 
strained to uniformly spaced values of /b using a force 
constant of 20 a.u. The regions of /b = 0.0 — 0.05 and 
/b = 0.095 — 1.00 are each sampled with 51 uniformly 
spaced windows using a harmonic restraint force constant 
of 3000 a.u. All sampling trajectories are performed us- 
ing PIMD with the masses of the classical particles set 
to nig = M = 0.01 a.u; as before, the altered masses in 
the PIMD sampling trajectories allow for larger timesteps 
while having no effect on calculated ensemble averages. 
Each sampling trajectory is run for at least 12.09 ps, the 
PIMD timestep is 2.42 x lO"** fs, and thermostatting is 
performed by resampling velocities from the MB distri- 
bution every 2.42 fs. The FE profiles are constructed 
from the sampling trajectories using WHAM. 

For each value of e, the transmission coef- 
ficient in Model SBl is calculated from 2400 
RPMD trajectories released from the dividing sur- 
face and evolved for 121 fs with the timestep of 
1.21 X lO^"' fs. The position of the dividing surface is 
fl = (0.5, 0.385, 0.2345, 0.014, 0.014, 0.014) for the 
e-cases (I, II, ..., VI). In Model SB2, 1600 RPMD 
trajectories are released at each value of e; each 
trajectory is evolved for 121 fs using a timestep of 
2.42 X lO"** fs; and the dividing surface is located at 
fl = (0.5, 0.65, 0.75, 0.986, 0.986, 0.986) for e-cases 
(I, II, ..., VI). Initial configurations for the released 
RPMD trajectories are sampled every 14.5 fs from eight 
long, independent PIMD sampling trajectories that are 
constrained to the dividing surface. The velocities of 
the PIMD sampling trajectories are resampled every 
48.4 fs from the MB distribution. The dividing surface 
constraint is implemented using the RATTLE algorithm. 

We note that RPMD results can be affected by cou- 
pling of fictitious internal ring-polymer modes to physi- 
cal frequencies in the system. -- We thus performed test 
calculations of the ET rate in these and similar sys- 
tems using pa rtially adiabatic centroid molecular dynam- 
ics (PACMD):22lI]The PACMD calculations revealed no 
significant changes from the RPMD results, confirming 
that this issue does not impact our conclusions. 



For the calculation of Marcus theory rates, the reorga- 
nization energy and the thermodynamic driving force for 
each value of epsilon are obtained analytically from the 
diabatic states for the donor and acceptor, Vii(s) and 
V22(s). For Model SBl, we obtain a solvent reorganiza- 
tion energy of A = 68.9 kcal/mol, and for Model SB2, we 
obtain A = 17.0 kcal/mol. The values of |Vi2| for both 
system-bath models is given in Table |IV[ 



3. Semiclassical Instanton Theory 

For the SB models, contributions from the linearly- 
coupled harmonic bath can be factorized and cancelled 
from the RHS of Eq. [l2] yielding expressions that de- 
pend only on the electron ring-polymer coordinates and 
the single classical solvent coordinate, s. Calculation of 
fcsci then consists of (i) determination of saddle-point 
configurations for the classical action, A ({g^"-*}; s) , on a 
numerical grid in the solvent coordinate s, (ii) evaluation 
of the steepest-descent approximation for I„(s) at each 
point on the solvent grid, and (in) integration over the 
solvent coordinate in Eq. [19] via numerical quadrature. 
The reactant partition function, Qr, was similarly ob- 
tained by evaluating /„(s) via steepest-decent expansion 
around the minimum-action configuration in the reactant 
basin. All calculations were performed using n — 2048 
beads for the electron ring polymer. 

For Model SBl, the grid in the solvent coordinate s 
consists of 200 uniformly spaced points in the range of 
—4 to 4 a.u.; for Model SB2, this grid consists of 150 uni- 
formly spaced points in the range —3 to 3 a.u. At each 
value of s, the saddle-point configuration on the surface 
A ({q^"'}; s) corresponds to the maximum along the path 
of minimum action that connects the reactant and prod- 
uct basins. This path of minimum action is obtained 
using the string method, ^^ with the path discretized 
into L = 1000 equidistant slices and with minimiza- 
tion performed using Euler integration and a timestep of 
2.4 X 10^^ fs. Initial convergence of the path is achieved 
when this minimization results in a change of less than 
5.3 X 10~* A in each degree of freedom. The path is then 
iteratively refined in the vicinity of the saddle point: a 
20-slice sub-section of the path about the saddle point 
is extracted, the number of slices used to describe the 
path is doubled, and the sub-section of the path is re- 
minimized with its endpoints fixed. Iterative refinement 
of the path is complete when the slice of maximum ac 
tion (i.e., the saddle point configuration) satisfies Eq. 16 
to within 10~^ a.u. 



4. QUAPI 

The QUAPI calculation for Model SB2 requires con- 
struction of the short-time system propagator followed 



10 



by two independent Monte Carlo (MC) simulations to 



evaluate the flux-flux correlation function in Eq. 29 



The complex-time propagator in Eq. [33] is calculated 
using eigenvalues and eigenfunctions obtained from a 2D 
discrete variable representation (DVR) grid calculation22l 
in the solvent coordinate, s, and the electronic state 
variable, a. The DVR Hamiltonian is diagonalized on 
a grid of 40 uniformly spaced points over a range of 
—4 to +4 a.u. in s and a = 1,2. The number of eigen- 
values and eigenvectors used in these calculations (Mq in 
Eq. [33]) ranges from 30 to 50 for the values of e considered 
in this study. 

The flux-flux correlation function in Eq.[29]is obtained 
from standard path-integral Monte Carlo (PIMC) sam- 
pling performed on the 2D DVR grid. In a first PIMC 
simulation, the correlation function is obtained using 

CpFit) = Dp{ sgn{Re[/i(s,cr;ie) -/2(s,cr;tc) 

+ /3(s, a: Q - his, a- U)]} )p(,_^^,^) , (60) 

where importance sampling is performed using the dis- 
tribution 

p(s, a; tc) = Abs {Re [/is((t; t^) ~ his, cr; tc) 

+ 4(s,cr;ic) -4(s,cr;te)]}, (61) 



and the function /^(s, cr; tc) is defined in Eq. 31 Conver- 
gence is achieved with 10^ MC steps. The normalization 
constant, Dp, is obtained from a second, independent 
PIMC simulation, using 



Do^Da 



Pis,cr;h) 

A(s,(T;ic) 



(62) 



A(s,(T;tc) 



Here, importance sampling is performed on the distribu- 
tion 



A(s,(T;ic)= Y\ \{'^k,Sk\e 
k=N'+3 

AA+l 

xn 

fc=2 



HsAt'/h 



Wk-l,Sk-l) 



^fc,.Sfe|e-*«^^*=/^|(Tfc_i,Sfe_i; 



(63) 



where cti = 2, ctaa+i = 2, a^f+2 = 1, and a2Af+2 = 1- 
Convergence is achieved with 10^ MC steps, and the 
normalization constant Dj^ is obtained by direct matrix 
multiplication. A maximum of A/" = 4 path beads are re- 
quired to converge the flux-flux correlation function over 
a timescale of 25 fs; no significant changes are observed 
between calculations performed using Af — A and Af — 8. 
The reactant partition function is obtained from a sin- 
gle PIMC calculation using the expression 



QK = TT[e-^"Vi], 



(64) 



where Pi = |1)(1| is the projection operator for the reac- 
tant electronic state. 

The QUAPI calculations for case IV are performed us- 
ing a larger value for the coupling between the solvent co- 
ordinate and the bath modes, -q/Mwc = 30. This change 



TABLE VI. ET reaction rates for the atomistic representa- 
tion, obtained using RPMD and Marcus theoryF] 

Case — AG° logfcMx log/cRPMD 



I 


0.0 


-2.2(1) 


-1.7(2) 


II 


22(1) 


4.6(3) 


4.6(4) 


III 


43(1) 


8.7(2) 


8.4(2) 


IV 


63(1) 


10.40(4) 


10.17(9) 


V 


84(1) 


10.0(1) 


11.21(4) 


VI 


124(3) 


2.9(8) 


11.48(8) 


VII 


138(2) 


-1.8(9) 


11.80(7) 



ET rates are given in s^^ , and the —AG" are given in 
kcal/mol. The numbers in parentheses denote the statistical 
uncertainty in the last reported digit. 



leads to lower-amplitude oscillations in the flux-flux cor- 
relation function and improved numerical convergence 
of the ET rate calculation. Other features of the flux- 
flux correlation function, including the timescales for the 
real-time oscillations and the decorrelation time, are un- 
changed. The invariance of these features suggests that 
the parameters used in the current study correspond to 
the regime in which the ET reaction rate is indepen- 
dent of the solvent-bath coupling.®'^ ^^ RPMD rate calcu- 
lations performed using different values for -q/Muc also 
support this conclusion. 



V. RESULTS 

A. Atomistic Simulations 

The atomistic representation for ET (Fig. IT]) is inves- 
tigated using direct RPMD simulations and the Marcus 
rate theory. For each case of the thermodynamic driv- 
ing force. Fig. l2Fa) presents FE profiles for the reactant 
and product diabatic electronic states as a function of 
the solvent collective variable, Af/(Q) (Eq. I59|. The 
FE profiles are obtained by reducing the corresponding 
2D surfaces, -F(/b, AC/), where the reactant and product 
diabats are associated with ring-polymer configurations 
for which /b > 0.995 and /b < 0.005, respectively. The 
results in Fig. l2ja) are graphically identical to those ob- 
tained using the tight-binding approximation, and the 
FE profiles exhibit the anticipated parabolic form, al- 
though no assumptions rega rding the linear response of 
the solvent have been madcl^SEll These data, in combi- 
nation with the calculated tunnel splitting for the trans- 
ferring electron, are used to calculate the Marcus rates 
in Table Ell 

Fig. [2lb) presents the corresponding FE profiles as 



a function of the bead-count coordinate, /b (Eq. 58 1. 
These profiles are used in the statistical component of 
the RPMD rate calculation (Eqs. Is] and [4]). As is seen 
from the inset, all of the profiles behave similarly in the 
vicinity of /b ~ 1. The steep rise in the FE profile be- 
tween 0.980 and 0.999 is associated with the formation 
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FIG. 2. (a) FE profiles, F(A(7), for tfie reactant (colored, at 
riglit) and product (left) diabatic electronic states as a func- 
tion of the solvent collective variable in the atomistic rep- 
resentation. The various cases of the thermodynamic driving 
force for the ET reaction are labeled; see Tables [TTT] and [Vl] for 
details. For each case, the FE profiles are vertically shifted to 
align the minima of the product basin, (b) The correspond- 
ing FE profiles as a function of the bead-count coordinate, 
F{fh). In the main panel, the profiles are vertically shifted to 
align the product basin; in the inset, the profiles are vertically 
shifted to align the reactant basin, (c) The corresponding 
RPMD transmission coefficients for the ET reaction, K{t). In 
panels (b) and (c), the curves retain the same color scheme 
introduced in panel (a). 



of "kink-pair" configurations, in which the ring polymer 
spans both redox sites, ^^ ^^ ^^ a typical kink-pair config- 
uration is illustrated in Fig. IT|a). 

The dynamical component of the RPMD rate calcu- 
lation (Eq. 11) is obtained from the long-time platearP^ 
of the RPMD transmission coefficient shown in Fig. W[c) . 
Each transmission coefficient is calculated with respect to 
a dividing surface at a fixed value of /b, as is described 
in Sec. |IV A| Plateau values in the range of 0.1-0.4 indi- 
cate modest recrossing of the RPMD trajectories through 
these surfaces. For cases in which the thermodynamic 
driving force corresponds to ET in the normal and aciti- 
vationless regimes. Fig. W[c) illustrates that the RPMD 
trajectories commit to the reactant or product basins 
within 10-20 fs, the timescale for local solvent motion be- 
tween librational rebounds. At thermodynamic driving 
forces corresponding to the inverted regime, the transmis- 
sion coefficient plateaus on faster timescales than those 
involving the rigid solvent molecules. 

Fig. Isja) presents a direct comparison of the RPMD 
and Marcus theory rates throughout the normal and acti- 
vationless regime for ET in the atomistic representation. 
The RPMD rates, which are also reported in Table VI 



quantitatively agree with the Marcus theory results over 
12 orders of magnitude in the ET reaction rate. Unlike 
the Marcus rates, which are based on a TST description 
for the reaction, the calculated RPMD rates are indepen- 
dent of any a priori assumptions about the ET reaction 
mechanism. 

Figs.lslb) and (c) illustrate the ET reaction mechanism 
that is predicted from the RPMD simulations. Represen- 
tative RPMD trajectories are projected onto the {z, AU) 
plane, where z is the component of the ring-polymer cen- 
troid that lies along the axis of the metal ions in the sys- 
tem; also shown are FE profiles for the system in these 
collective variables. For symmetric ET (Case I), Fig.lS^b) 
reveals that the RPMD trajectories involve three distinct 
steps that will be familiar from the Marcus rate theory: 
(i) solvent fluctuation to a configuration for which the 
reactant and product diabats are nearly degenerate (in- 
dicated by the dashed line) , (ii) formation of a kink-pair 
in the ring-polymer configuration and rapid transfer of 
the electron from one redox site to the other, and (Hi) 
relaxation of the solvent coordinate in the product basin 
following the ET event. For ET approaching the activa- 
tionless regime (Case IV), Fig. ISlc) shows the latter two 
steps in the mechanism remain, but only a small initial 
solvent fluctuation is needed to reach solvent configura- 
tions for which the electronic diabats are degenerate. 

To understand the connection between RPMD and the 
Marcus theory rate expression, we note that Eq. 37 in- 
cludes two key terms - an Arrhenius-type contribution 
that is associated with free energy of solvent reorganiza- 
tion to bring reactant and product diabats into degener- 
acy and a prefactor that depends on the coupling between 
the diabatic states. RPMD captures the solvent reorgani- 
zation energetics because the path- integr al-based method 
preserves exact quantum statistics.'22E3 The RPMD rate 
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FIG. 3. (a) ET reaction rates for the atomistic representa- 
tion in the normal and activationless regimes, computed using 
RPMD (red) and Marcus theory (black). The various cases 
for the thermodynamic driving force are labeled, (b) Rep- 
resentative trajectories (red) from the ensemble of reactive 
RPMD trajectories for symmetric ET (Case I). The trajec- 
tories are plotted as a function of the ring-polymer centroid, 
z, and the solvent collective variable, At/. The FE profile in 
these collective variables is also presented, with contour lines 
indicating FE increments of 10 kcal/mol. (c) Representative 
RPMD trajectories for activationless ET (Case IV) and the 
corresponding FE profile. The white arrows in panels (b) 
and (c) indicate the solvent reorganization mechanism for ET 
that is anticipated in the Marcus rate theory, and the dashed 
lines indicate values of A(7 at which the reactant and product 
diabats cross in Fig. [2ra) . 



also correctly accounts for the tunneling contribution to 
the ET reaction rate, which can likewise be attributed to 
the path-integral basis of the method; the tunnel split- 
ting for the electron between degenerate redox sites is 
analytically related to the reversible work for form ing a 
kink-pair in the ring-polymer configuration.ESESMl Given 
that the ensemble of reactive RPMD trajectories exhibit 
the dual rare events of solvent reorganization and kink- 
pair formation, and given that the FE barriers associated 
with these two steps are analytically related to the key 
terms in the Marcus rate expression, it is reasonable that 
Fig. IsFa) finds good agreement between RMPD and Mar- 
cus theory. The RPMD method succeeds in the normal 
and activationless regimes because it captures the correct 
physics of the ET reaction. 

Fig. |4] demonstrates that the success of the RPMD 
method does not extend into the inverted regime for ET, 
with both the RPMD rates and reaction mechanism devi- 
ating from the predictions of Marcus theory. In Fig. l4[a) , 
the RPMD rates are seen to be only weakly dependent 
on the increasing driving force, rather than exhibiting 
the characteristic turnover in this inverted regime. The 
RPMD trajectories also deviate from the reaction mech- 
anism that is assumed in the Marcus TST, as is seen 
in Fig. Qb). The reactive trajectories exhibit kink-pair 
formation directly from solvent configurations that are 
characteristic of the reactant basin; the expected solvent 
reorganization to configurations for which the electronic 
diabats are degenerate (indicated by the dashed line in 
the figure) is not observed. 

To further explore the successes and failures of RPMD 
in these various regimes for ET, we compare the method 
with semiclassical instanton theory and exact quantum 
dynamics in the following section. 



B. System-Bath Simulations 

In this section, we employ system-bath representations 
for ET to allow for the comparison of RPMD with other 
simulation techniques, including semiclassical instanton 
and exact quantum dynamics methods. 

Fig. Isja) and Table Vll present a comparison of the 
RPMD and Marcus rates for Model SBl, which is pa- 



rameterized to match the energy-scales for the atom- 
istic representation (Sec. IllB). As before, the RPMD 



method reproduces the Marcus rates throughout the nor- 
mal and activationless regimes, while failing to predict 
the turnover of the ET rate in the inverted regime. Anal- 
ysis of the RPMD reactive trajectories in this system re- 
veals mechanisms that are entirely analogous to those 
observed in Figs. [sFb), [slc), and Wb) for the atomistic 
system. Specifically, for the normal and activationless 
regimes, the RPMD trajectories exhibit solvent reorga- 
nization to configurations for which the electronic diabats 
are degenerate, followed by rapid transfer of the electron 
between redox sites; and for the inverted regime, RPMD 
predicts ET without prior solvent reorganization. These 
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FIG. 4. (a) ET reaction rates for the atomistic representation 
in the inverted regime, computed using RPMD (red) and Mar- 
cus theory (black). The various cases for the thermodynamic 
driving force are labeled, (b) Representative trajectories (red) 
from the ensemble of reactive RPMD trajectories for inverted 
ET (Case VI). The trajectories are plotted as a function of 
the ring-polymer centroid, z, and the solvent collective vari- 
able, A[/. The FE profile in these collective variables is also 
presented, with contour lines indicating FE increments of 10 
kcal/mol. The white arrows indicate the solvent reorganiza- 
tion mechanism for ET that is anticipated in the Marcus rate 
theory, and the dashed line indicates the value of A(7 at which 
the reactant and product diabats cross in Fig. [2[a). 



data confirm that Model SBl exhibits the same essential 
physics as the atomistic representation. 



ET rates from the steepest-descent SCI theory (Eq. 19 ) 
are also included in Fig. IsFa) and Table VII Through- 
out the full range of thermodynamic driving forces, the 
instanton method tracks the RPMD results, including 
deviation from the Marcus predictions in the inverted 
regime. As is shown in Table |VIH a- correct ion of the 
RMPD rates (Eq. 22 assuming Kq ~ 1) further im- 
proves their agreement with the SCI rates. These results 



FIG. 5. (a) ET reaction rates for Model SBl, computed using 
RPMD (red), Marcus theory (black), and SCI theory (blue), 
(b) FE profiles, F(A[/), for the reactant (colored, at right) 
and product (left) diabatic electronic states as a function of 
the solvent coordinate, s. The various cases of the thermody- 
namic driving force for the ET reaction are labeled; see Table 
IIVI for details. The arrow indicates the value of the solvent 
coordinate that maximizes Iji(s), which corresponds to the 
dominant contribution to the SCI rate in Eq. |19| 



TABLE VII. ET reaction rates for Model SBl, obtained using 
RPMD, Marcus theory, and SCI theory^ 



]ase 


-AG« 


log fcMT 


log fcRPMD 


log fcsci 


log afcRPMD 


I 


0.0 


-6.0 


-6.55(4) 


-7.9 


-7.3 


II 


18.5 


-0.2 


-0.33(3) 


-1.9 


-0.9 


III 


36.9 


3.7 


3.52(8) 


2.4 


2.8 


IV 


73.9 


6.3 


6.19(5) 


5.6 


5.4 


V 


110.9 


1.6 


6.44(1) 


5.8 


5.8 


VI 


148.0 


-10.4 


6.69(3) 


5.9 


5.9 



ET rates are given in s ^, and the —AG" are given in kcal/mol. 
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TABLE VIII. ET reaction rates for a ID asymmetric double 
well, obtained using SCI theory^ 

e AS log fcsci 



0.0 
0.05 
0.10 
0.20 
0.30 
0.40 



0.0 
0.02940 
0.05884 
0.11776 
0.17676 
0.23584 



-11.0 
-10.9 
-10.8 
-10.6 
-10.3 
-10.3 



The Golden Rule for the symmetric case yields logfc = —11.55. 
AE is the difference between the two lowest eigenenergies for 
the system. All quantities reported in atomic units. 



underscore that the failure of RPMD does not arise from 
a breakdown in its formal connection with SCI theory)^ 
instead, the comparison suggests that both RPMD and 
the SCI theory share the same underlying flaw in the 
inverted regime .123 

The mechanistic predictions from SCI theory also show 
similarities with the RPMD results. Fig. Isjb) presents 
the Marcus parabolas for the electronic diabats of Model 
SBl as a function of the solvent coordinate, s. Also 
shown are the solvent configurations that correspond to 
the SCI predictions for the ET transition state. For each 
value of the thermodynamic driving force, the arrow in 
the figure indicates the solvent configuration that maxi- 
mizes I„(s), which corresponds to the largest contribu- 
tion to the rate in Eq. [T9j For the normal and acti- 
vationless regimes, SCI theory correctly predicts an ET 
transition state at the crossing of the electronic diabats. 
However, in the inverted regime, the SCI transition state 
is instead located at the minimum of the reactant basin. 
These mechanistic results from SCI theory are consistent 
with the observed pathways for the RPMD trajectories, 
which suggests that in the inverted regime, both RPMD 
and SCI theory overestimate the degree of ET from sol- 
vent configurations in the reactant basin. 

To further illustrate this issue, we present SCI rate 
calculations for deep tunneling in a ID asymmetric dou- 
ble well. Table |VIII| presents ET reaction rates calcu- 
lated on the potential energy surface Uo-m{q) (Eq. 49 1, 
with parameters from Model SBl. Although this is a 



non-dissipative ID system, the SCI rate is still well- 
defined, and it is reported as a function of the poten- 
tial energy asymmetry. The rates plateau to a finite 
value with increasing asymmetry, which is consistent 
with rates for dee p tu nneling between a bound state 
and a continuum.'^SHiool However, this behavior is qualita- 
tively incorrect for tunneling rates between bound states, 
which should vanish for n on-d egenerate states in accord 
with Fermi's Golden Rule.ES^We conclude that SCI the- 
ory, as well as the closely related RPMD method, signif- 
icantly overestimate the tunneling probability between 
asymmetric bound states, leading to an incorrect ET 
mechanism and overestimation of the reaction rate in the 
inverted regime. 
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FIG. 6. The ET rates for Model SBl corresponding to a 
Marcus-like mechanism (black) and the "direct" mechanism 
in Eq. 66 (red). SCI rates (blue) correspond to the kinetically 



favorable mechanism in all regimes. See text for details. 



The results for the simple double-well system can be 
used to deduce a more general argument for the applica- 
bility of the RPMD and SCI calculations in ET problems. 
Table |VIII combined with the condition of detailed bal- 
ance for the thermal reaction rate, indicates that the SCI 
rate for transfer in an asymmetric double-well system is 
approximately 



k 



27r 



T/i2p min(l,e" 



0AE\ 



(65) 



For the Marcus-type ET mechanism in which electron 
tunneling is gated via solvent reorganization that sym- 



metrizes the double- well system, Eq. 65 leads to the TST 
rate in Eq. [37] However, for an unphysical "direct" ET 
mechanism in which electron tunneling proceeds from sol- 
vent configurations in the reactant basin (i.e., without 
prior solvent reorganization), Eq. 65 leads to the follow- 



ing TST expression for the ET rate. 



Wct = ^|T^i2|^(^) mm 



in I 



l,e 



-P{x+ag'') 



(66) 



Fig. [6] presents the ET reaction rates for Model SBl, 
assuming either the Marcus-type mechanism (Eq. 37 
black) or the direct mechanism (Eq. l66l red). Also 
plotted are the rates calculated using SCI theory 
(Eq. 19 blue). Throughout the normal and activationless 



regimes, the rate for the Marcus-type mechanism domi- 
nates; in the inverted regime, the rate for the direct mech- 
anism dominates; and the results from SCI theory closely 
track the larger of these two rates. It is clear that SCI 
theory (as well as RPMD) features a competition between 
the correct, Marcus- type mechanism for ET and the un- 
physical, direct mechanism for ET, and the prevailing 
mechanism is that which is predicted to be faster. This 
analysis is fully consistent with the earlier discussions of 
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function (Eq. 25) contributes to the turnover in the ET 



2 4 6 8 10 12 14 16 

t/fs 



reaction rate in the inverted regime. The RPMD approx- 
imation to the real-time dynamics of the system, which 
is not expected to capture coherent quantum effects J22HiJ 
does not fully enforce the quantization of electronic dy- 
namics and leads to the observed inaccuracies in the in- 
verted regime. Approximate quantum dynamical meth- 
ods that explicitly enforce electronic quantization by us- 
ing either a discrete electronic state basis or by exactly 
mapping to a continuous electronic basis are thus ex- 
pected to provide a better starting point for describing 
state-t o-st atc electronic dynamics and ET in the inverted 
regime.Q2Ii4 i6 19 104 pm-ther investigation of this point is 
in progress. 



FIG. 7. Normalized flux-flux autocorrelation functions CFp(t) 
for Model SB2, calculated using exact quantum dynamics for 
Cases I (black), II (blue). III (purple) and IV (red). 



TABLE IX. ET reaction rates for Model SB2, obtained us- 
ing RPMD, Marcus theory, SCI theory, and exact quantum 
dynamics F] 



Case 


-AG° 


log fcMT 


log feRPMD 


log fcsci 


logfcq 


I 


0.0 


6.7 


6.05(3) 


5.1 


6.7(1) 


II 


5.3 


8.3 


7.73(5) 


6.6 


8.5(1) 


III 


8.8 


9.1 


8.54(3) 


7.4 


9.0(3) 


IV 


17.6 


9.8 


9.27(2) 


8.6 


10.8(9) 


V 


26.5 


8.9 


9.40(3) 


8.8 


- 


VI 


35.3 


6.3 


9.52(2) 


8.8 


— 



ET rates are 
kcal/mol. 



given in s ^, and the —AG" are given in 



the RPMD trajectories (Figs. [3)d, |3|:, and ^) and the 
SCI transition state configurations (Fig. l5pj for ET in 
the various regimes. Furthermore, this analysis provides 
a general basis for expecting the SCI and RPMD meth- 
ods to accurately describe ET rates in the normal and 
activationless regimes, and for expecting these methods 
to significantly overestimate the ET rate in the inverted 
regime. 



Table IX presents ET rates for Model SB2, including 
results obtained using the QUAPI exact quantum dy- 
namics method. Comparison of the RPMD, Marcus the- 
ory, and SCI theory rates for ET in Table |IX| confirms 
that Model SB2 exhibits all of the previously discussed 
trends for these approximate methods. Fig. [7] presents 
the flux-flux correlation functions used to obtain the ex- 
act quantum rates for Model SB2. 

The results in Fig. [7] emphasize the role of elec- 
tronic state quantization in the ET reaction dynamics. 
At larger thermodynamic driving forces, the correlation 
functions become increasingly oscillatory, with a reso- 
nance frequency that matches the electronic state e nergy 
gap between the ET reactant and product .E^SEnH Inte- 
gration over this increasingly oscillatory time correlation 



VI. CONCLUSIONS 

The current paper demonstrates the applicability of 
RPMD for the direct simulation of ET reaction dynamics 
in complex systems. Using both atomistic and system- 
bath representations for ET in a polar solvent, we com- 
pare RPMD results with those obtained using Marcus 
theory, semiclassical instanton theory, and exact quan- 
tum dynamics. Throughout the normal and activation- 
less regimes for ET, RPMD correctly predicts the ET 
reaction mechanism and quantitatively describes the ET 
reaction rate over 12 orders of magnitude, without in- 
voking any prior mechanistic or transition state assump- 
tions. Analysis of the RPMD trajectories reveals that 
the accuracy of the method lies in its exact description of 
statistical fluctuations, with regard to both solvent reor- 
ganization and the formation of kink-pair configurations 
during the electron tunneling event. However, for ET in 
the inverted regime, both RPMD and SCI theory fail to 
predict the turnover in the ET reaction rate with increas- 
ing thermodynamic driving force. In this regime, both 
methods overestimate the probability of electronic tun- 
neling from solvent configurations in the reactant basin, 
leading to an overestimation of the corresponding reac- 
tion rates. Exact quantum dynamics calculations illus- 
trate that the limitations of the RPMD method in the 
inverted regime arise from the inadequate quantization of 
the real-time electronic-state dynamics; analogous break- 
downs of the method have been identified in other appli- 
cations to strongly coherent quantum swtems, includ- 
ing low-dimensional quantum oscillators^S and electron- 
scattering in dilute fiuids.**^ 

We conclude by emphasizing that the normal and acti- 
vationless regimes encompass the vast m ajori ty ET reac- 
tions in biological and synthetic systems.l^S3 -pj^g results 
presented here thus constitute a significant success for 
the RPMD method, demonstrating that it allows for the 
robust, direct simulation of thermally activated ET in 
systems with over 1000 atoms, leading to the quantita- 
tive prediction of ET reaction rates and and the poten- 
tial discovery and characterization of ET reaction mecha- 
nisms in complex systems. A comparable demonstration 
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using other approximate real-time quantum simulation 
methods has not, to our knowledge, been previously re- 
ported. Having established both the applicability and 
limitations of RPMD for ET reactions dynamics, this 
work provides the foundation for future studies of ET 
and proton-coupled ET reactions in enzymes and other 
condensed-phase systems. 
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Appendix A: System-Bath Potential Energy Parameters 



TABLE 
electron- 
SB2. 



Xn. Parameters for the left Coulombic well in the 
ion potential energy function of Eq. [49] for Model 



Case 


fflL 


6l 


cl 


^L" 


out 


I 


0.157480 


1.596286 


1.609286 


-3.050000 


-7.086414 


II 


0.157507 


1.596671 


1.612042 


-3.050000 


-7.087150 


III 


0.157525 


1.596927 


1.613880 


-3.050000 


-7.087642 


IV 


0.157569 


1.597568 


1.618474 


-3.050000 


-7.088871 


V 


0.157613 


1.598208 


1.623065 


-3.050000 


-7.090102 


VI 


0.157657 


1.598848 


1.627656 


-3.050000 


-7.091336 



Unless otherwise noted, parameters are given in atomic units 



TABLE XIII. Parameters for the right Coulombic well in the 
electron-ion potential energy function of Eq. |49] for Model 
SB2. 



Case 


Or 


6h, 


CR 


r'S 


'^R 


I 


0.157480 


-1.596286 


1.609286 


3.050000 


7.086414 


II 


0.156666 


-1.587919 


1.598482 


3.050000 


7.085675 


III 


0.156124 


-1.582341 


1.591280 


3.050000 


7.085179 


IV 


0.154767 


-1.568396 


1.573273 


3.050000 


7.083924 


V 


0.153410 


-1.554450 


1.555265 


3.050000 


7.082650 


VI 


0.152053 


-1.540504 


1.537256 


3.050000 


7.081357 



^ Unless otherwise noted, parameters are given in atomic units 



TABLE X. Parameters for the left Coulombic well in the 
electron-ion potential energy function of Eq. [49] for Model 
SBl.[f] 



Case 


flL 


feL 


Cl 


^L 


rT 


I 


0.164567 


2.002721 


3.683127 


-4.062912 


-8.106702 


II 


0.164520 


2.001856 


3.675494 


-4.062912 


-8.104948 


III 


0.164472 


2.000989 


3.667859 


-4.062912 


-8.103197 


IV 


0.164377 


1.999251 


3.652576 


-4.062912 


-8.099709 


V 


0.164280 


1.997506 


3.637280 


-4.062912 


-8.096237 


VI 


0.164183 


1.995754 


3.621971 


-4.062912 


-8.092782 



Unless otherwise noted, parameters are given in atomic units 



TABLE 
electron- 
SBl.[f] 



XL Parameters for the right Coulombic well in the 
ion potential energy function of Eq. [49] for Model 



Case 


flR 


&R 


cr 


r'S 


out 
'^R 


I 


0.164567 


-2.002721 


3.683127 


4.062912 


8.106702 


II 


0.167357 


-2.036963 


3.752141 


4.062912 


8.108432 


III 


0.170147 


-2.071204 


3.821152 


4.062912 


8.110110 


IV 


0.175726 


-2.139680 


3.959165 


4.062912 


8.113319 


V 


0.181304 


-2.208150 


4.097166 


4.062912 


8.116346 


VI 


0.186882 


-2.276615 


4.235157 


4.062912 


8.119207 



Unless otherwise noted, parameters are given in atomic units 



Appendix B: Transformation to Diabatic Basis 

The QUAPf method is implemented in an electronic 
diabatic state representation of the ET reaction. In this 
appendix, we describe the procedure for transforming the 
potential energy function for Model SB2 from a position 
basis for the electron (Eq. 



48) to a diabatic basis where 



the reactant and product electronic states are maximally 
localized on the donor and acceptor metal atoms. 

We begin by calculating the two lowest adiabatic elec- 
tronic eigenstates {ipo{q;s) and tpiiq^s)) and eigenen- 
ergies (-Eo(s) and Ei{s)) of the system Hamiltonian at 
fixed values values of the solvent coordinate in the range 
—8 ap < s < 8 ap. For each value of s, the system Hamil- 
tonian is diagonalized on a uniform DVR grid of 1024 
electron positions in the range —25 ao < <Z < 25 oq. 

For each value of s, reactant and product electronic 
wavefunctions in the diabatic basis are obtained via ro- 
tation of the two lowest-energy adiabatic wavefunctions, 
using 

(f>R{q; s) = cos{9s)'ipo{q; s) - sin(6's)V'i((7; s) (Bl) 

and 

(j)p{q; s) = sm{0s)'4'o{q; s) + cos{es)i^i{q; s), (B2) 

where 



1 , / Sio + Sqi 
- arctan ^ 



(B3) 



and Sfj,i, = J_ il>fj,{q;s)*ipi^{q',s)dq. This choice of 
the rotation angle, 9g, maximizes J_ |0r('?;'S)| dq, the 



17 



TABLE XIV. The diagonal elements of the diabatic potential 
matrix Vii(s) and V22(s) in Eqs. 55 and 56 for Model SB2. 



Case ai x 10^ 


6i X 10^ 


I 4.7722 


1.1308 


II 4.7722 


1.1308 


III 4.7722 


1.1308 


IV 4.7720 


1.1307 



Ci 02 X lO'' 62 X 10-^ C2 

-2.1576 4.7722 1.1308 -2.1576 

-2.1477 4.7722 1.1308 -2.1561 

-2.1411 4.7721 1.1308 -2.1551 

-2.1245 4.7720 1.1308 -2.1526 



probability that the reactant diabatic state is positioned 
on the donor ion. Maximization of the probability that 
the product diabatic state is positioned on the acceptor 
ion yields an identical choice for Og. 

The corresponding potential energy matrix elements in 
the diabatic basis (Eq. 23 1 are thus 



Vii(s) = Eo{s) cos2 Os + Ei{s) sin^ 0„ (B4) 

V22is)^Eo{s) sin^ 61, + Ei{s) cos^ Os , (B5) 

V2i{s) ^ Vi2{s) = {Eo{s)-Ei{s)) cos 6s smOs. (B6) 

The diagonal elements are found to be parabolic func- 
tions of s, and the off-diagonal element are found to be 
nearly constant with respect to s. We fit Vii(s) and 
^22(5) to second-order polynomials functions (Eqs. |55] 
and 56 ) and employ a constant value for V12 that corre- 



sponds to the s ~ result. The polynomial expansion co- 



efficients for Vii(s) and V22(s) are provided in Table XIV 



and the constant value for V12 is provided in Table IV 
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